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NUMERICAL  SIMULATION  OF  SOLAR  CORONAL  MAGNETIC  FIELDS 


I.  INTRODUCTION 

It  is  widely  believed  that  the  solax  magnetic  field  is  the  driver  for  essentially  all  man¬ 
ifestations  of  coronal  activity.  Both  the  structure  and  the  dynamics  of  coronal  plasma  are 
determined  by  the  magnetic  field.  The  basic  idea  is  that  photospheric  footpoint  motions 
stress  the  coronal  field  lines,  thereby  determining  the  coronal  topology  and  producing  mag¬ 
netic  free  energy  that  can  be  converted  into  plasma  kinetic  and  thermal  energy.  Hence,  the 
response  of  the  coronal  magnetoplasma  to  footpoint  motions  is  of  fundamental  importance 
for  understanding  solar  activity  and  many  other  astrophysical  phenomena. 

In  recent  years  there  has  been  a  great  deal  of  work  on  this  topic  of  coronal  response  to 
photospheric  motions.  Due  to  the  complexity  of  the  relevant  equations,  most  of  the  effort 
has  been  on  numerical  simulations.  We  present  in  this  paper  the  results  of  a  numerical 
simulation  that  applies  to  the  situation  of  twisting  a  flux  tube  in  a  magnetic  arcade.  The 
calculation  is  a  dynamic,  fully-three-dimensional  MHD  simulation  in  a  geometry  appropri¬ 
ate  to  the  corona.  To  our  knowledge  this  is  the  first  one  of  its  kind.  The  model  and  the 
numerical  methods  are  described  in  the  next  Section. 
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The  results  of  our  simulation  are  pertinent  to  several  outstanding  questions  in  solar 
coronal  studies.  One  question  concerns  the  topology  of  prominence  magnetic  fields.  It  is 
commonly  believed  that  the  cool,  dense,  H0  emitting  material  in  high  quiescent  promi¬ 
nences  must  be  supported  by  the  magnetic  field  (e.g.  Tandberg-Hanssen  1974).  In  that 
case  the  magnetic  field  lines  must  form  a  “hammock”  shape  in  the  corona  so  that  the 
cool  material  cam  rest  at  the  bottom  of  the  hammock.  Several  hammock- type  models  have 
been  proposed  (e.g.  Kippenhahn  and  Schliiter  1957,  Kuperus  and  Raadu  1974).  The  key 
question  is,  how  does  the  field  obtain  such  a  geometry?  A  single  current-free  arcade  has 
field  lines  that  are  concave-down  in  the  corona.  Perhaps  the  simplest  and  most  appealing 
explanation  is  that  photospheric  motions  twist  the  field  lines  of  an  initially  current-free 
flux  tube  into  the  appropriate  shape.  Parker  (1979)  has  shown  that  the  twist  will  tend  to 
propagate  to  the  weakest  part  of  the  field,  near  the  apexes  of  the  flux  tube.  If  sufficient 
twist  is  introduced,  (of  order  one  complete  rotation  or  more),  and  if  the  twist  is  sufficiently 
localised  at  the  apex,  then  some  of  the  field  lines  must  acquire  a  concave-up  geometry. 
Analytic  models  for  such  a  twisted  flux  tube  have  been  constructed  by  severed  authors 
(e.g.  Demoulin  and  Priest  1989).  However,  the  analytic  models  consider  only  an  equilib¬ 
rium  flux  tube  in  isolation  and  only  for  linear  departures  from  a  2-dimension?  1  equilibrium. 
They  do  not  consider  the  effect  of  a  realistic  arcade  structure  for  the  field,  the  dynamic 
formation  of  the  twists  and  the  effect  of  instabilities.  Our  simulation  allows  us  to  study 
all  these  topological  effects. 
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The  role  of  ideal  instabilities  is  an  impoitm  issue  in  models  for  coronal  activity. 
Since  many  observed  phenomena  such  as  prominence  eruption,  sprays  and  coronal  mass 
ejections  exhibit  mass  motions  with  velocities  of  the  order  of  the  coronal  Alfven  speed,  it 
seems  likely  that  they  involve  an  ideal  instability  such  as  a  kink.  The  basic  idea  is  that  if 
the  coronal  field  is  stressed  sufficiently  by  photospheric  motions,  it  becomes  unstable  and 
erupts  outward  carrying  coronal  and  chromospheric  plasma  with  it  (e.g.  Sturrock  1989). 
Two  recent  3-d  simulations  seem  to  support  this  scenario  (Strauss  and  Otani  1988;  Mikic 
1990).  These  authors  found  that  continuous  twisting  of  a  flux  tube  eventually  resulted 
in  a  kink  instability.  However,  their  simulations  were  performed  for  a  geometry  in  which 
the  field  is  confined  between  two  vertical  plates  and  the  system  has  cylindrical  symmetry 
(before  kinking).  These  constraints  permit  the  buildup  of  arbitrarily  large  stresses.  In  the 
true  coronal  geometry  the  field  is  free  to  expand  upward  and  there  is  no  symmetry  to  the 
field.  Hence  the  coronal  magnetic  field  can  relax  gradually  in  response  to  photospheric 
motions.  It  is  not  clear  if  kinking  will  occur  in  this  case  as  well,  hence,  more  simulations 
are  needed  to  determine  whether  ideal  instabilities  occur  in  more  realistic  coronal  magnetic 
fields. 

Another  issue  is  the  role  of  magnetic  reconnection.  A  number  of  authors  (e.g.  Pneu- 
man  1983,  Moore  1988,  Sturrock  1989,  Van  Ballegooijen  and  Martens  1989,  etc)  have 
invoked  reconnection  in  a  sheared  arcade  as  a  mechanism  for  both  the  formation  and 
eruption  of  prominences.  The  concept  underlying  their  models  is  that  reconnection  can 
redistribute  the  shear  in  an  arcade  such  that  part  of  the  arcade  flux  becomes  less  sheared 
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while  another  part  becomes  progressively  more  sheared.  In  their  model  reconnection  has 
the  effect  of  concentrating  the  shear  in  a  certain  flux  region.  It  is  this  highly  sheared  flux 
that  eventually  erupts,  and  leads  to  flares  (Moore  1988).  This  is  a  very  promising  model. 
Simulations  on  a  reduced  two-dimensional  model,  coupled  with  the  assumption  that  the 
reconnection  did  rearrange  the  shear  as  postulated,  do  show  prominence  eruption  (Vein 
Ballegooijen  and  Martens  1989).  The  full  three-dimensional  model  has  yet  to  be  verified 
by  rigorous  calculations.  However,  there  are  serious  doubts  as  to  whether  true  reconnec¬ 
tion  (as  opposed  to  simple  diffusion)  can  ever  take  place  in  a  coronal  arcade.  Linear  theory 
(e.g.,  Velli  and  Hood  1986)  indicates  that  resistive  tearing  is  suppressed  by  line-tying  in  a 
field  that  does  not  have  a  reversal  of  the  unsheaxed  initial  field,  as  is  the  case  for  a  dipole 
arcade.  In  the  general  non-linear  regime,  Greene  (1988)  has  presented  arguments  for  the 
absence  of  reconnection  in  a  field  with  the  topology  of  a  single  arcade.  Hence,  we  believe 
that  the  role  of  reconnection  in  a  coronal  arcade  is  an  open  question  which  detailed,  3-d 
simulations  can  help  to  resolve. 

A  related  issue  concerns  the  formation  of  DC  electric  current  sheets  in  closed  magnetic 
configurations  in  the  solar  corona.  Current  sheets  are  very  important  for  models  of  coronal 
heating.  In  order  to  account  for  the  observed  energy  losses,  it  is  widely  believed  that  the 
electric  currents  in  the  corona  must  attain  a  highly  filamented  structure  as  in  a  current 
sheet.  Parker  (1972,  1987)  proposes  that  current  sheets  form  as  a  natural  consequence 
of  footpoint  motions.  He  argues  that  a  general  3-dimensional  magnetic  field  such  as  that 
in  the  corona  can  have  no  well-behaved  equilibrium  state,  and  that  magnetic  singularities 
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corresponding  to  current  sheets  must  form  to  allow  the  magnetic  field  to  relax  to  a  lower 
energy  state.  Several  authors  have  numerically  tested  this  hypothesis  using  3-dimensional 
force-free-field  codes  (Van  Ballegooijen  1988,  Mikic  et  al,  1989).  They  concluded  that 
their  calculations  did  not  support  Parker’s  hypothesis.  These  calculations,  however,  were 
performed  for  an  initially  uniform  field  confined  between  two  plates.  Our  arcade  geometry 
is  more  appropriate  for  the  corona,  and  we  include  the  dynamics,  so  that  we  can  examine 
different  aspects  of  the  current-sheet  formation  question  them  the  equilibria  codes. 

In  the  next  section  we  describe  the  exact  model  that  we  use  for  the  simulation  and 
our  numerical  code.  The  final  section  discusses  the  results  and  our  conclusions. 

II.  THE  SIMULATION  MODEL 


a)  The  Numerical  Code 

We  begin  with  the  dissipative,  incompressible  magnetohydrodynamic  (MHD)  equa¬ 
tions,  written  in  a  dimensionless  rotation  form: 


dv  „  1 

—  =  v  x  u  -  VII  +  j  x  B  +  77 V2v, 
at  M 


(la) 


and, 


V  •  v  =  0, 

d  A  _  _ ,  1  _  _ 

-5—  =  vxB  +  V$  -  —  V  x  V  x  A. 
at  b 

V  •  B  =  0, 


(16) 

(lc) 


(Id) 
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where  v(x,  t)  =  ( u,v,w )  =  flow  velocity,  u>(x,t)  =  V  x  v  =  vorticity,  A(x,<)  =  magnetic 
vector  potential,  B(x,t)  =  V  x  A  =  magnetic  field,  j(x,  t)  =  V  x  B  =  electric  current 
density,  II(x,t)  =  mechanical  pressure  +  kinetic  energy  density,  $(x,  t)  =  magnetic  scalar 
potential,  S  =  Ltmdquist  number,  and  M  =  viscous  Lundquist  number.  Lundquist  num¬ 
bers  are  Reynolds-like  numbers,  with  the  Alfven  speed  used  for  the  characteristic  velocity. 
We  choose  a  gauge  with  $  =  0.  In  this  representation,  B  is  measured  in  terms  of  a  char- 
acteristic  field  strength,  e.g.,  Bq  =  (F-1(|  B  |  ,  where  the  brackets  <>  indicate  an 

integral  over  the  system  and  V  denotes  the  volume  of  integration.  The  velocities  are  mea¬ 
sured  in  units  of  the  Alfven  speed,  Ca  =  Bo/(4np)i  where  the  mass  density,  p,  is  assumed 
to  be  constant  and  uniform.  For  a  characteristic  distance,  Lq,  lhe  time  is  measured  in 
units  of  the  Alfven  transit  time,  Lq/Ca ■  The  Lundquist  number  S  =  CaL0/t?,  where  r) 
is  the  magnetic  diffusivity.  The  viscous  Lundquist  number  M  =  CaLq/v,  where  v  is  the 
kinematic  viscosity. 

For  the  simulation  reported  here,  rj  —  .0025  and  v  =  .2  so  that  5  =  500  and  M  =  5. 
The  low  viscous  Lundquist  number  was  chosen  to  damp  out  fluid  motions.  In  addition, 
it  helps  to  dissipate  small-scale  magnetic  structure  via  the  Alfven  effect,  allowing  us  to 
perform  fully  3D  numerical  simulations  at  a  relatively  large  magnetic  Lundquist  number. 
However,  such  small  values  for  the  Lundquist  numbers  might  be  thought  to  preclude  mag¬ 
netic  reconnection  a  priori.  To  establish  the  possibility  of  reconnection  at  these  Lundquist 
numbers,  we  have  examined  several  linear  and  nonlinear  problems  where  we  have  previ¬ 
ously  observed  magnetic  reconnection.  One  such  case  is  the  linear  magnetic  reconnection 
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problem,  which  reduces  to  the  MHD  version  of  the  Orr-Sommerfeld  equation  (Dahlburg 
et  al.,  1983,  1986).  For  the  mean  magnetic  field  Bx(z)  =  tanh(z),  we  find  that  linearly 
unstable  modes  do  exist  at  these  Lundquist  numbers.  In  terms  of  the  wavenumber  parallel 
to  the  mean  magnetic  field,  a,  and  the  growth  rate,  u>,-,  we  find,  for  example;  (a,u;,)  = 
(0.5,  0.01),  (1.0,  0.028057),  (2.0,  0.03878),  (3.0,  0.02255),  and  (4.0,  0.0096).  This  implies 
that  linear  instabilities  are  at  least  possible  at  these  Lundquist  numbers.  We  have  also 
considered  a  highly  nonlinear  initial  value  reconnection  problem,  the  Orszag-Tang  vortex, 
at  these  Lundquist  numbers  (Orszag  and  Tang  1979,  Dahlburg  and  Picone  1989).  Re¬ 
connection  occurs  within  less  than  one  Alfven  transit  time,  with  the  attendent  formation 
of  electric  current  sheets  and  vortex  quadrupoles,  in  a  manner  similar  to  previously  re¬ 
ported  numerical  simulations  of  this  configuration.  Hence  we  also  conclude  that  nonlinear 
magnetic  reconnection  can  occur  at  these  Lundquist  numbers.  We  have  also  run  parts 
of  the  magnetic  arcade  simulation  reported  in  the  present  paper  with  a  smaller  viscosity, 
M  =  30,  and  have  found  no  significant  qualitative  change  in  the  results.  Let  us  emphasize 
that  these  Lundquist  numbers,  although  large  by  the  standards  of  numerical  simulations, 
are  still  many  orders  of  magnitude  smaller  than  the  expected  solar  values;  hence,  the 
simulations  are  bound  to  be  much  more  diffusive  than  the  real  corona. 

Equations  1  a^e  solved  in  a  three-dimensional  Cartesian  geometry.  We  employ  the 
following  numerical  algorithm  to  discretize  the  equations:  A  Fourier  pseudospectral  dis¬ 
cretization  is  used  in  x  and  y,  the  periodic  directions.  Chebyshev  collocation  is  imple¬ 
mented  in  the  z  direction,  with  the  mechanical  pressure  evaluated  on  a  grid  staggered  in  z. 
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The  temporal  discretization  is  a  low-storage,  third-order,  Runge-Kutta  -  Crank-Nicolson 
scheme,  with  a  backward-Euler  pressure  correction  at  each  time-level  (Zang  and  Hussaini 
1986).  For  a  grid  of  322  Fourier  modes  and  32  Chebyshev  modes  the  code  requires  about 
2  seconds  per  time  step  on  the  NRL  CRAY  X-MP,  and  approximately  .8  million  words  of 
memory.  More  details  on  the  algorithm  are  given  in  the  appendix. 

b)  The  Initial  and  Boundary  Conditions 

Since  the  code  requires  periodicity  in  the  x  and  y  directions  we  pick  as  the  initial  field 
a  periodic  array  of  current-free  arcades.  The  initial  vector  potential  is  given  by: 

Ax  =  Az  =  0,  (2a) 

Ay  —  cosfcxexp(— k(z  +  1))  ,  (26) 

and  the  computational  domain  consists  of  0  <  x  <  2n,0  <  y  <  2n,—  1  <  z  <  1.  This  is 
the  same  potential  field  studied  in  2.5  dimensions  by  Mikic  et  al.,  1988  who  considered 
the  evolution  of  this  field  under  the  influence  of  photospheric  shear  flows.  The  initial  field 
has  translational  symmetry  along  the  y  direction;  however,  this  symmetry  is  broken  im¬ 
mediately  by  the  flow  imposed  at  the  bottom  boundary  of  the  computational  box,  which 
is  supposed  to  correspond  to  the  photosphere.  Unlike  previous  calculations  (Mikic  et  al., 
1989,  Van  Ballegooijen  1988)  which  focused  on  equilibrium  states,  we  concentrate  on  coro¬ 
nal  dynamics  and  do  not  attempt  to  obtain  an  equilibrium  or  steady-state.  Our  simulation 
consists  of  two  distinct  periods,  a  “driving”  part  during  which  free  energy  is  built  up  in 
the  coronal  magnetic  field  by  footpoint  motions,  and  a  “decay”  part  during  which  the  field 
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relaxes  back  to  a  current-free  state.  The  field  is  first  twisted  with  a  prescribed  footpoint 
motion  for  10  characteristic  time  scales  r,  where  r  is  the  average  Alfven  crossing  time  for 
the  whole  system.  Then  the  footpoint  positions  are  fixed,  and  the  system  is  allowed  to 
evolve  for  another  60r.  During  the  driving  period,  t  <  lOr,  the  horizontal  velocities  at  the 
bottom  wall  z  =  —  1  are  assumed  to  be, 


dip  dip 

V,  =  —  and  a,  = 


where  ip  is  given  by, 


ip  =  C  sin2  2kx  sin3  ky , 


(3) 

(4) 


and  C  is  a  constant.  Of  course  during  the  decay  period,  t  >  lOr,  vx  =  vy  =  0.  Note 
that  the  form  (3)  for  the  footpoint  flow  corresponds  to  an  incompressible  twisting  motion. 
The  bottom  wall  is  also  assumed  to  be  impenetrable,  so  that  vz  =  0  throughout  the 
simulation.  The  boundary  conditions  at  the  other  walls  are  straightforward.  At  the  sides 
of  the  computational  box  we  must  assume  periodic  boundary  conditions.  The  top  wall, 
z  =  1,  is  assumed  to  be  a  rigid,  impenetrable  surface,  consequently  v  =  0  there. 


The  only  other  quantity  that  requires  boundary  conditions  is  the  vector  potential  A. 
Again  periodic  conditions  are  imposed  at  the  sides.  The  top  and  bottom  conditions  on  A 
are  selected  so  that  the  evolution  at  the  boundary  corresponds  to  ideal  flow.  However,  this 
determines  only  the  horizontal  components  of  A, 

~dt  ~  VyBz'  (5<l) 
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(5b) 


dAy 

dt 


Vz  B  zi 


where  we  have  used  the  fact  that  vz  vanishes  on  the  boundary.  In  principle,  no  boundary 
conditions  need  be  imposed  on  Az ,  but  since  the  numerical  code  requires  them,  we  simply 
solve  the  partial  differential  equation  for  Az  at  the  boundary  ( cf .  Mikic  et  al.,  1988). 


dAz 

dt 


=  (v  x  B),  - 


(5c) 


We  have  also  tried  extrapolating  its  value  from  the  interior,  but  the  results  axe  not  sensitive 
to  this  assumption. 


Due  to  the  combination  of  the  periodicity  restriction  and  the  particular  problem  under 
consideration,  there  is  an  unavoidable  redundancy  in  the  calculations.  There  are  two 
arcades  in  the  computational  box,  one  full  arcade  with  its  neutral  line  at  the  center  of 
the  box,  x  =  27r/6,  and  one  half  arcade  on  either  side.  Throughout  the  simulation  we 
found  that  both  arcades  behave  identically,  (as  they  should);  hence,  we  will  only  show  the 
central  arcade  in  the  following  figures.  This  redundancy  will  not  be  present  in  calculations 
with  more  realistic  magnetic  fields  and  photospheric  flows.  For  example,  the  photospheric 
flow  pattern  is  perhaps  best  represented  as  a  random,  time  evolving  flow-field.  Such  a 
photospheric  flow  pattern  would  break  the  symmetry  implied  by  equations  2  and  4. 

The  value  of  C  in  equation  (4)  for  the  twisting  motion  is  chosen  so  that  the  maximum 
velocity  at  the  photosphere  is  approximately  1/10  of  the  coronal  Alfven  speed.  Of  course, 
this  ratio  for  the  driving  speed  to  the  characteristic  velocity  is  much  greater  than  in  the 
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solar  case.  Typical  photospheric  velocities  are  of  order  1  km/s  and  coronal  Alfven  velocities 
~  1, 000  km/s  so  that  the  actual  ratio  should  be  of  order  10-3.  However,  it  is  not  possible 
to  use  such  a  low  ratio  in  the  numerical  simulations  because  diffusion  is  so  much  stronger 
in  the  simulation  model  than  on  the  real  Sun.  If  the  driving  speed  is  too  low,  then  the 
effects  of  the  footpoint  motions  will  diffuse  away  before  any  twist  can  be  built  up  in  the 
coronal  magnetic  field;  the  field  is  able  to  “slip”  back  faster  than  it  is  driven.  Also,  since 
the  time  step  in  the  code  is  set  by  the  characteristic  velocity,  it  is  simply  too  expensive 
computationally  to  use  a  very  low  ratio.  A  value  of  ~  1/10  is  about  the  smallest  that  can 
be  used  with  our  present  code.  This  implies  that  our  results  during  the  driving  period  are 
likely  to  be  too  dynamic  compared  to  the  solar  corona;  however,  the  simulation  should  be 
accurate  for  the  decay  phase. 

In  Figure  la  the  stream  lines  of  the  photospheric  driving  flow,  equivalent  to  con¬ 
stant  contours  of  ip,  are  shown,  and  in  Figure  lb,  contours  of  the  velocity  magnitude, 
| v |  =  [(dip/dx )2  -F  (dip/dy)2]  3.  The  flow  pattern  has  been  chosen  so  that  the  velocity 
vanishes  at  the  boundary  between  adjacent  arcades  and  at  the  neutral  line  of  each  ar¬ 
cade.  This  ensures  that  each  arcade  retains  its  identity;  there  is  no  mixing  of  flux  between 
different  arcades  or  between  opposite  polarity  regions.  From  Figure  1  we  note  that  the 
velocities  are  quite  concentrated  and  are  oppositely  directed  on  either  side  of  the  neutral 
line  so  that  they  correspond  to  twisting  a  fairly  localized  flux  region.  Again  due  to  the  pe¬ 
riodicity  requirement  there  is  a  redundancy  in  the  velocity  pattern,  consequently  we  have 
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two  regions  of  twisted  flux.  They  behaved  identically  throughout  the  simulation,  which 
lends  confidence  in  the  numerical  code. 

The  footpoint  displacements  that  result  from  the  photospheric  velocities  are  shown  in 
Figure  2,  where  we  plot  the  contours  of  constant  Bz  on  the  bottom  wall  at  the  beginning 
and  the  end  of  the  driving  period,  t  =  lOr.  Note  that  from  equation  (2)  the  initial 
photospheric  flux  Bz  =  —ksinkx  so  that  the  initial  contours  are  straight  lines  parallel  to 
the  y  axis  (Figure  2a).  From  the  ideal  induction  equation  one  finds  that  Bz  is  convected 
with  the  flow  for  an  incompressible  horizontal  motion, 

dBz 

-^+vPrVB,  =  0,  (6) 

where  v^,/,  is  the  horizontal  velocity  at  the  photosphere.  Hence,  by  following  contours  of 
Bz  one  can  track  the  footpoint  positions. 

There  are  several  important  features  of  Figure  2.  It  clearly  demonstrates  the  point 
that  for  ideal  motions  the  topology  of  the  coronal  field  can  be  determined  from  the  bound¬ 
ary  conditions  at  the  photosphere.  The  fact  that  the  coronal  field  has  been  twisted  is 
obvious  from  the  photospheric  contours.  Even  if  the  footpoint  motions  were  much  more 
complicated,  involving  winding  or  braiding  patterns  and  compressible  flow,  the  topology 
would  be  recoverable  from  the  contours  of  the  Euler  potentials  for  B  at  the  photosphere 
(Antiochos  1987,1990).  It  is  also  evident  from  Figure  2  that  due  to  the  footpoint  motions, 
the  Bz  contours  have  been  “stretched  out”,  and  become  more  concentrated.  We  conclude 
that  the  gradients  in  the  footpoint  displacements  increase  with  time  and,  therefore,  the 
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scales  of  the  magnetic  structures  in  the  corona  must  decrease  with  continued  twisting  (e.g. 
Van  Ballegooijen  1988,  Mikic  et  al.,  1989).  Note  that  this  occurs  even  if  the  scales  of  the 
photospheric  velocities  remain  constant.  This  result  implies  that  a  photospheric  velocity 
with  large  scales  such  as  that  given  by  (4),  will  eventually  produce  structure  on  the  photo¬ 
sphere,  (and  consequently  in  the  corona),  with  scales  below  the  numerical  resolution.  The 
actual  magnitude  of  the  twisting  is  evident  in  Figure  2.  Each  side  of  the  arcade  has  been 
twisted  through  a  maximum  angle  ~  180°,  so  that  the  total  twist  is  approximately  one 
full  rotation,  which  is  about  the  maximum  twist  that  we  cam  resolve  with  a  323  grid.  In 
order  to  resolve  a  more  complicated  motion  such  as  a  braiding  pattern  which  requires  at 
least  three  rotations  (Antiochos  1990),  we  would  need  a  much  larger  grid. 

III.  RESULTS  OF  THE  SIMULATION 

The  evolution  of  the  magnetic  field  lines  during  the  simulation  are  shown  in  Figures 
3a  -  3d  (side  view),  and  4a  -  4d  (top  view).  We  plot  two  sets  of  field  lines  at  four  times 
during  the  run,  t  =  0,  lOr,  20r,  and  70r.  These  field  lines  were  chosen  by  the  position  of 
one  of  their  footpoints  on  the  photosphere.  We  selected  two  rows  of  17  points  lying  on  the 
photosphere,  and  then  traced  the  field  line  emanating  from  these  points.  The  two  rows 
are  parallel  to  the  y  axis  and  pass  approximately  through  the  center  of  the  twisted  region 
of  the  arcade.  In  the  Figures,  these  points  are  to  the  left  of  the  photospheric  neutral  line 
which  lies  at  the  center  of  the  grid,  (x,  =  in).  Field  lines  with  footpoints  on  the  inner 
row  are  shown  as  dotted  lines,  those  on  the  outer  are  solid. 
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Figures  3a  and  4a  show  the  initial  current-free  magnetic  field  (equation  2).  The 
translational  symmetry  along  the  y  direction  is  apparent.  The  field  clearly  has  a  very 
simple  topology  with  all  of  the  field  lines  from  the  inner  row  of  footpoints  lying  well  under 
the  lines  from  the  outer  row.  Figures  3b  and  4b  shows  the  field  at  the  end  of  the  driving 
phase,  t  =  lOr,  when  it  has  the  maximum  twist.  It  should  be  noted  that  these  are  not  the 
same  field  lines  as  in  3a  and  4a.  The  field  lines  of  Figures  3a  and  4a  have  moved  from  their 
original  positions  so  that  their  leftside  footpoints  no  longer  lie  on  our  two  straight  rows  (see 
figure  2).  However,  there  is  nothing  special  about  the  particular  field  lines  in  3a  and  4a, 
and  we  find  it  easier  to  see  the  effect  of  the  motions  by  keeping  a  fixed,  two-straight-rows 
pattern  for  the  leftside  footpoints  in  all  the  figures. 

Comparing  3a, 4a  and  3b, 4b  it  is  evident  that  the  field  topology  becomes  much  more 
complex  as  a  result  of  the  twisting.  The  field  lines  now  wrap  around  each  other  so  that 
there  is  no  clear  distinction  between  an  outer  and  an  inner  set  on  the  right  hand  side  of  the 
photospheric  neutral  line.  The  most  striking  feature  of  the  twisted  field  lines  is  that  they  are 
all  greatly  expanded  upward  from  the  current-free  positions.  This  result  is  expected  from 
studies  of  force-free  fields  (e.g.,  Aly  1984;  Yang,  Sturrock  and  Antiochos  1986).  The  field 
acts  to  minimize  any  stress  due  to  footpoint  motions  by  expanding  outward  in  the  corona. 
Note  that  this  effect  is  not  incorporated  in  models  with  an  initially  uniform  field  confined 
between  two  horizontal  plates  (e.g.,  Parker  1972,  Mikic  et  ai,  1989).  For  a  given  footpoint 
motion,  such  models  will  yield  fields  with  finer  current  structures  than  in  a  more  realistic 
open  coronal  geometry.  In  fact,  our  results  also  underestimate  the  outward  expansion. 
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The  boundary  conditions  that  we  assume  at  the  top  and  bottom  of  an  impenetrable  wall, 
and  the  restriction  to  incompressible  flows,  equation  (lb),  all  act  to  suppress  an  outward 
expansion. 

The  expansion  plays  a  crucial  role  in  determining  the  field’s  structure  and  dynamics. 
We  note  from  Figures  3b  and  4b  that  all  the  field  lines  are  concave  down  in  the  corona. 
There  are  no  “dips”  suitable  for  supporting  cool  prominence  material.  If  the  field  lines 
had  not  expanded  upward,  we  believe  that  dips  would  have  formed.  The  reason  for  the 
lack  of  dips  is  that  outward  expansion  tends  to  decrease  the  pitch  of  field  lines  as  they 
twist  around  each  other.  It  seems  likely  that  more  twisting  will  eventually  produce  dips, 
on  the  other  hand,  Aly  (1984)  has  argued  that  continual  stressing  of  the  field  should  lead 
to  an  open  arcade.  If  so  then  even  infinite  twisting  will  not  produce  dips  in  the  corona. 
Unfortunately,  our  simulation  is  not  capable  of  accurately  treating  more  than  about  one 
rotation  so  that  we  cannot  address  this  interesting  question  at  present. 

The  topology  of  the  field  lines  at  later  times  is  shown  in  Figures  3c  -  3d  and  4c  - 
4d.  Since  the  footpoints  are  held  fixed  during  this  phase  of  the  evolution,  then  these 
field  lines  would  correspond  exactly  to  those  in  3b, 4b  if  the  evolution  were  ideal.  It  is 
evident  from  the  Figures  that  the  evolution  is  far  from  ideal.  The  right-side  footpoints  are 
clearly  changing  with  time  even  though  we  allow  no  flow  on  the  boundary.  The  reason  for 
this  is  resistive  diffusion,  which  permits  the  field  lines  to  “slip”  back  to  their  current-free 
positions.  It  is  important  to  emphasize  that  this  slippage  is  not  a  boundary  effect;  it  is 
occurring  throughout  the  corona.  On  the  contrary,  at  the  boundary  the  code  is  able  to 
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ensure  an  ideal,  perfectly  rigid  plate;  but  in  the  interior  volume  the  code  must  use  a  finite 
resistivity.  Note  that  this  slipping  is  not  unphysical;  it  also  occurs  in  the  reeil  corona,  but 
at  a  much  slower  rate. 

We  find  that  the  evolution  throughout  the  decay  phase  is  one  of  simple  diffusion. 
It  is  evident  from  Figures  3  and  4  that  all  the  field  lines  tend  to  untwist.  There  is  no 
evidence  for  any  type  of  instability,  ideal  or  resistive.  This  can  be  seen  from  plots  of  the 
global  quantities.  We  show  in  Figure  5  the  evolution  during  the  whole  simulation  of  the 
logarithm  of  the  kinetic  energy  (Ev  =  l/2(jv|2}),  the  magnetic  energy  ( Eb  =  1/2(|B|2)), 
and  the  mean  square  electric  current  (J  =  l/2(|j|2)).  The  forcing  at  the  wall  causes  Ev  to 
build  rapidly  to  a  saturation  level.  The  kinetic  energy  decreases  rapidly,  due  to  the  high 
viscosity,  at  the  end  of  the  driving  phase  ( t  =  10)and  remains  small  thereafter.  Growth  of 
Ev,  typically  seen  when  jets  form  as  a  result  of  reconnection  (e.g.  Dahlburg,  Dahlburg,and 
Marislca  1988),  is  not  observed.  The  magnetic  energy,  Eb,  also  grows  quickly  during  the 
forcing  phase  to  about  1.5  times  its  initial  value.  After  the  forcing  is  cut  off  at  t  =  10,  Eb 
exhibits  a  smooth  monotonic  decay  to  a  level  below  its  initial  value.  This  reflects  the  fact 
that  the  magnetic  field  relaxes  to  a  different  equilibrium,  consistent  with  the  new  footpoint 
positions  exhibited  in  Figure  2b.  In  the  same  way,  J  increases  during  the  forcing  phase, 
as  the  forcing  creates  smaller  spatial  scale  structures  in  the  magnetic  field.  During  the 
damping  phase,  J  decays  to  near  zero,  indicating  the  relaxation  of  the  magnetic  field  to  a 
nearly  potential  state. 
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The  nature  of  the  magnetic  field  evolution  can  be  described  in  terms  of  an  effective 


magnetic  Lundquist  number,  given  by 


S(t)  = 


2([il2) 

£<|B|2). 


(7) 


This  is  also  shown  in  Figure  5.  We  obtain  S(t )  by  considering  the  magnetic  induction 
equation  in  the  limit  of  negligible  velocity,  taking  the  scalar  product  of  this  equation 
with  B,  and  integrating  the  resulting  magnetic  energy  density  equation  over  the  system. 
Solving  for  the  Lundquist  number  gives  the  desired  result.  Figure  5d  shows  that  S(t) 
returns  rapidly  to  the  predicted  Ohmic  rate,  viz.,  S(t)  =  500,  when  the  forcing  stops  at 
T  =  lOr.  The  “jiggling”  in  the  data  is  due  to  small  errors  in  computing  the  denominator 
of  S(t),  which  is  postprocessed.  Since  the  magnetic  energy  fails  to  decay  at  a  higher  rate, 
magnetic  reconnection  must  be  absent. 


Another  important  result  of  our  simulation  is  that  it  shows  no  evidence  for  current 
sheet  formation.  Hence,  our  results  are  consistant  with  the  studies  of  force-free  equilibria 
(Sakurai  1979;  Van  Ballegooijen  1988;  Mikic  et  al.,  1989).  The  lack  of  current  sheet 
formation  can  be  seen  from  the  evolution  of  J  shown  in  Figure  5.  Not  only  does  the 
integrated  current  exhibit  a  monotonic  decay,  but  we  have  also  examined  in  detail  the 
spectrum  of  the  current  density,  and  have  found  that  the  spectrum  decays  as  would  be 
expected  for  pure  resistive  diffusion. 
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IV.  DISCUSSION 


The  results  presented  in  this  paper  are  of  a  negative  character,  :.e.,,  we  do  not  see 
magnetic  reconnection,  kinking,  or  the  formation  of  concave-up  magnetic  field  lines  suitable 
for  prominence  formation.  Hence,  it  appears  that,  contrary  to  many  models,  photospheric 
twisting  of  a  single  arcade  does  not  lead  to  the  type  of  processes  required  to  explain  solar 
activity.  The  key  question  is  whether  this  conclusion  is  valid  in  general,  or  whether  it  is  due 
solely  to  the  limitations  of  this  particular  simulation.  The  numerical  simulation  itself  has 
three  major  shortcomings.  First,  it  is  too  diffusive,  especially  when  compared  to  the  real 
corona.  There  is  little  that  can  be  done  about  this  problem.  Larger  Lundquist  numbers  can 
be  obtained  by  going  to  larger  grids,  but  there  is  no  hope  at  present  for  obtaining  numbers 
significantly  larger  than  ~  1000.  Second,  the  numerical  resolution  is  barely  adequate.  At 
the  time  of  maximum  twist,  10r,  the  current  spectrum  shows  a  decrease  of  only  4  orders  of 
magnitude  over  the  wavenumber  range.  Ideally,  we  would  like  to  have  at  least  a  6  orders  of 
magnitude  decrease.  Since  we  have  used  only  a  323  grid,  and  much  larger  grids  are  possible 
on  present  computers,  we  see  no  difficulty  in  obtaining  substantially  better  resolution  in 
future  simulations.  A  third  problem  is  the  incompressibility  assumption.  This  is  not  valid 
for  the  solar  corona;  but,  we  believe  it  has  not  had  a  significant  effect  on  the  results.  If 
anything,  a  compressible  simulation  should  exhibit  an  even  larger  upward  expansion  and 
less  of  a  tendency  to  instability. 
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There  are  also  limitations  with  the  assumed  form  for  the  photospheric  motions.  We 
believe  that  the  magnitude  of  the  twisting  is  sufficient.  It  is  evident  from  Figure  5b  that 
the  magnetic  free  energy  is  50%  the  energy  of  the  potential  field.  Also  it  can  be  seen 
from  Figures  3b  and  4b  that  the  field  lines  are  strongly  twisted.  Coronal  loops  rarely 
appear  to  be  so  distorted.  Although  the  magnitude  of  the  motions  is  acceptable,  their 
size  scale  is  limited.  Due  to  the  difficulty  with  obtaining  adequate  numerical  resolution, 
we  have  considered  only  a  large-scale  twisting  in  which  the  width  of  the  twisting  motion 
is  approximately  equal  to  the  width  of  the  arcade.  This  situation  would  correspond  to  a 
large-scale  shearing  of  a  coronal  arcade,  and  is  the  easiest  to  resolve  numerically.  However, 
the  photosphere  also  contains  small-scale  motions  such  as  granular  motions,  for  example.  It 
may  be  that  the  small-scale  motions  play  a  crucial  role  in  processes  such  as  the  formation 
of  prominence  dips.  Simulations  with  larger  grids  are  needed  in  order  to  address  this 
possibility. 

Even  though  our  simulation  is  clearly  limited,  we  believe  that  at  least  some  of  the 
results  will  hold  for  general  footpoint  motions  and  Lundquist  numbers.  In  particular, 
we  expect  that  for  a  dipolar  magnetic  arcade,  resistive  instabilities  will  not  occur  due 
to  the  effect  of  line-tying,  irrespective  of  the  footpoint  motions.  This  argues  against  the 
reconnection  scenarios  proposed  by  a  number  of  authors.  We  also  expect  that,  contrary 
to  the  non-equilibrium  model,  current  sheets  will  not  form.  Of  course  the  scale  of  the 
currents  in  the  corona  is  determined  by  the  scale  of  the  footpoint  displacements,  so  that 
if  the  displacements  have  a  very  fine  scale,  then  fine  scale  currents  will  result.  However, 
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our  simulation  and  all  the  equilibria  calculations  (Sakurai  1979;  Van  Ballegooijen  1988; 
Mikic  et  al.,  1989)  find  that  for  well-resolved  footpoint  displacements  smooth  current 
distributions  result. 

The  question  of  ideal  instabilities  and  the  formation  of  concave-up  field  lines  is  less 
clear.  We  find  that  both  are  suppressed  by  the  outward  expansion  of  the  field.  Such 
an  expansion  will  occur  for  any  form  of  footpoint  motions,  but  it  is  possible  that  strong, 
highly-localised  twisting  can  lead  to  the  formation  of  concave-up  field  lines  and  to  kink-like 
instabilities  in  small  regions  of  the  arcade.  If  so,  then  small  scale  motions  may  have  im¬ 
portant  consequences  for  prominence  formation  and  coronal  dynamics.  This  issue  requires 
further  study. 

The  general  conclusion  of  our  results  is  that  photospheric  stressing  of  a  magnetic  ar¬ 
cade  is  incapable  of  producing  phenomena  such  as  prominence  eruption  or  coronal  heating. 
Some  critical  ingredient  is  missing  from  this  model.  There  are  several  possibilities.  One  is 
that  emerging  flux  plays  the  dominant  role  in  activating  the  coronal  magnetic  field.  Sim¬ 
ulation  which  allow  for  flux  emergence  and  submergence  need  to  be  considered.  Another 
possibility  is  that  a  complex  magnetic  topology  is  required.  We  have  argued  (Antiochos 
1989,1990)  that  a  pure  dipolar  field  (i.e.,  one  due  to  a  single  positive  and  a  single  negative 
polarity  region  on  the  photosphere)  is  somewhat  pathological,  since  it  is  the  only  field  with 
a  well-behaved  connectivity  everywhere.  In  a  more  realistic  field  which  is  due  to  a  number 
of  different  polarity  regions  at  the  photosphere,  the  connectivity  must  be  discontinuous 


and  null  points  must  occur  in  the  corona.  In  this  case  current  sheets  and  reconnection  will 
readily  occur  as  a  result  of  footpoint  motions.  We  are  currently  pursuing  these  possibilities. 

In  summary,  we  believe  that  the  work  reported  here  constitutes  only  a  first  step  in 
examining  the  3D  dynamics  of  coronal  magnetic  fields.  Much  more  interesting  simulations 
await  to  be  done. 


APPENDIX 

COMPUTATIONAL  PROCEDURE 

We  modify  the  algorithm  of  Zang  and  Hussaini  (1986)  to  include  magnetic  field  terms. 
The  governing  equations  1  -  5  are  solved  within  the  computational  range  Q  by  means  of 
the  following  time-step  splitting: 


v*  =  v*  xw*+j*  xB*  +  -^-V2v  (in  Q) 

M 

At*  =  v*  x  B*  -  x  V  x  A  (in  ft) 

v*  =  g*  (on  dQ ) 

A*  =  h*  (on  dQ) 

v*(x, tn )  =  v(x,t„)  (in  ft) 
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A*(x,  tn)  =  A  (x,<„) 

(in  fl) 

v**  =  -VP** 

(in  fl) 

o 

II 

• 

» 

> 

O 

(in  Cl) 

v**  •  n  =  g  •  n 

(on  dCl ) 

v**  f  =  -VP**  •  f 

(on  dCl) 

v**(x,P)  =  v*(x,P) 

(in  Q) 

The  first  time-level,  represented  by  single-starred  variables,  is  a  convective,  dissipative 
step.  The  second  level,  represented  by  double-starred  variables,  is  a  pressure  step  which 
enforces  the  solenoidality  of  the  flow  field. 

A  staggered  grid  is  employed  using  the  following  collocation  points:  x,  =  2ni/Nz  for 
i  =  0,1  ,NZ  —  1,  j/j  =  2k j/Ny  for  j  =  0,1,  Ny  —  1,  r*  =  cos(nk/Nz )  for  k  =  0,1,  Nz, 
and,  Zfc+1/2  =  cos[7r(/:  +  1/2)/ ]  for  k  =  0, 1,NZ  —  1.  Velocities  and  vector  potentials  are 
defined  at  the  nodes  (x,,y;,2*)-  Pressures  are  defined  at  the  shifted  points  (xj,  yj,  Zk+i/2)- 
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PHOTOSPHERIC  FLOW  FIELD 
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v  field  given  in  equation  4,  with  C  =  1.  (a)  Stream  function 

■pe  neutral  line  in  the  center  of  the  flow  pattern,  running  along  the 

\Oi/dx)2  +  W/dy)2]2  . 


NORMAL  MAGNETIC  FIELD  AT  PHOTOSPHERE 
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Fig.  2  —  One  quadrant  of  the  normal  magnetic  field  in  the  photosphere.  The  same  quadrant  is  shown  in  Fig.  1.  (a) 
t  =  Or,  initial  conditions,  (b)  I  =  lOr,  conclusion  of  the  twisting  phase.  The  figure  shows  how  much  the  magnetic  field  is 
twisted  from  its  initial  configuration,  approximately  one  full  rotation.  After  t  =  lOr,  the  photospheric  B.  retains  this  form. 


Fig.  3  —  One  quadrant  of  magnetic  field  lines,  side  view.  The  same  quadrant  is  shown  as  in  Fig.  1.  (a)  t  =  Or,  initial 
conditions,  (b)  t  =  10r,  conclusion  of  the  twisting  phase,  (c)  t  =  20r,  partial  damping,  (d)  t  =  70r,  conclusion  of  the 
simulations. 


Fig.  4  —  One  quadrant  of  magnetic  field  lines,  top  view.  The  same  quadrant  is  shown  as  in  Fig.  1.  (a)  r  =  Or,  initial 
conditions,  (b)  t  =  lOr,  conclusion  of  the  twisting  phase,  (c)  t  =  20r,  partial  damping,  (d)  t  =  70r,  conclusion  of  the 
simulations. 


Fig.  5  -  Global  quantities  versus  time,  (a)  logarithm  of  kinetic  energy  (1/2  <  |  v  | :  >).  (b)  magnetic  energy 
( 1  /2  <  |  B| 2  > ) .  (c)  mean  square  electnc  current  ( 1  /2  <  |  j  | 2  > ) .  (d)  mean  magnetic  Lundquist  number 


